using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.QuadPack
{
    /// <summary>
    /// Integration rules for integrands having algebraico-logarithmic end point singularities.
    /// </summary>
    [Serializable]
    public static class Dqc25sClass
    {
        /// <summary>
        /// Integration rules for integrands having algebraico-logarithmic end point singularities.
        /// </summary>
        /// <param name="f">The f parameter.</param>
        /// <param name="a">The a parameter.</param>
        /// <param name="b">The b parameter.</param>
        /// <param name="bl">The bl parameter.</param>
        /// <param name="br">The br parameter.</param>
        /// <param name="alfa">The alfa parameter.</param>
        /// <param name="beta">The beta parameter.</param>
        /// <param name="ri">The ri parameter.</param>
        /// <param name="rj">The rj parameter.</param>
        /// <param name="rg">The rg parameter.</param>
        /// <param name="rh">The rh parameter.</param>
        /// <param name="abserr">The abserr parameter.</param>
        /// <param name="resasc">The resasc parameter.</param>
        /// <param name="wgtfunc">The wgtfunc parameter.</param>
        /// <param name="nev">The nev parameter.</param>
        /// <returns>For further information view the QuadPack documentation.</returns>
        public static double Dqc25s(HardRealFunction f, double a, double b, double bl, double br, double alfa,
                                    double beta, double[] ri, double[] rj, double[] rg, double[] rh, ref double abserr,
                                    ref double resasc, int wgtfunc, ref int nev)
        {
            double[] x = {
                             0.99144486137381041114, 0.96592582628906828675, 0.92387953251128675613,
                             0.86602540378443864676
                             , 0.79335334029123516458, 0.70710678118654752440, 0.60876142900872063942,
                             0.50000000000000000000, 0.38268343236508977173, 0.25881904510252076235,
                             0.13052619222005159155
                         };

            double centr;
            double dc;
            double factor;
            double fix;
            double hlgth;
            double resabs = 0;
            double res12;
            double res24;
            double u;
            double result;
            double[] cheb12 = new double[13];
            double[] cheb24 = new double[25];
            double[] fval = new double[25];
            int i;
            int isym;

            nev = 25;
            if ((bl == a) && ((alfa != 0.0) || (wgtfunc == 2) || (wgtfunc == 4)))
            {
                goto _10;
            }
            if ((br == b) && ((beta != 0.0) || (wgtfunc == 3) || (wgtfunc == 4)))
            {
                goto _140;
            }

            //  If a>bl and b<br, apply the 15-point Gauss-Kronrod scheme. 
            result = Dqk15wClass.G_K15W(f, DqwgtClass.Dqwgts, a, b, alfa, beta, wgtfunc, bl, br, ref abserr, ref resabs,
                                        ref resasc);
            nev = 15;
            goto _270;

            //  This part is only executed if a = bl.
            // *  Compute the Chebyshev series expansion of the following function:
            // *  f1 = (0.5*(b+b-br-a)-0.5*(br-a)*x)^beta*
            // *          f(0.5*(br-a)*x+0.5*(br+a))
            // 
            _10:
            hlgth = 0.5 * (br - bl);
            centr = 0.5 * (br + bl);
            fix = b - centr;
            fval[0] = 0.5 * f.SolveAt(hlgth + centr) * Math.Pow(fix - hlgth, beta);
            fval[12] = f.SolveAt(centr) * Math.Pow(fix, beta);
            fval[24] = 0.5 * f.SolveAt(centr - hlgth) * Math.Pow(fix + hlgth, beta);
            for (i = 1; i < 12; i++)
            {
                u = hlgth * x[i - 1];
                isym = 24 - i;
                fval[i] = f.SolveAt(u + centr) * Math.Pow(fix - u, beta);
                fval[isym] = f.SolveAt(centr - u) * Math.Pow(fix + u, beta);
            }
            factor = Math.Pow(hlgth, alfa + 1.0);
            result = 0.0;
            abserr = 0.0;
            res12 = 0.0;
            res24 = 0.0;
            if (wgtfunc > 2)
            {
                goto _70;
            }
            DqchebClass.Dqcheb(x, fval, cheb12, cheb24);

            //  wgtfunc = 1  (or 2) 
            for (i = 0; i < 13; i++)
            {
                res12 += (cheb12[i] * ri[i]);
                res24 += (cheb24[i] * ri[i]);
            }
            for (i = 13; i < 25; i++)
            {
                res24 += (cheb24[i] * ri[i]);
            }
            if (wgtfunc == 1)
            {
                goto _130;
            }

            //  wgtfunc = 2 
            dc = Math.Log(br - bl);
            result = res24 * dc;
            res12 = 0.0;
            res24 = 0.0;
            for (i = 0; i < 13; i++)
            {
                res12 += (cheb12[i] * rg[i]);
                res24 += (cheb24[i] * rg[i]);
            }
            for (i = 13; i < 25; i++)
            {
                res24 += (cheb24[i] * rg[i]);
            }
            goto _130;

            //  Compute the Chebyshev series expansion of the following function:
            // *      f4 = f1*log(0.5*(b+b-br-a)-0.5*(br-a)*x)
            // 
            _70:
            fval[0] *= Math.Log(fix - hlgth);
            fval[12] *= Math.Log(fix);
            fval[24] *= Math.Log(fix + hlgth);
            for (i = 1; i < 12; i++)
            {
                u = hlgth * x[i - 1];
                isym = 24 - i;
                fval[i] *= Math.Log(fix - u);
                fval[isym] *= Math.Log(fix + u);
            }
            DqchebClass.Dqcheb(x, fval, cheb12, cheb24);

            //  wgtfunc = 3  (or 4) 
            for (i = 0; i < 13; i++)
            {
                res12 += (cheb12[i] * ri[i]);
                res24 += (cheb24[i] * ri[i]);
            }
            for (i = 13; i < 25; i++)
            {
                res24 += (cheb24[i] * ri[i]);
            }
            if (wgtfunc == 3)
            {
                goto _130;
            }

            //  wgtfunc = 4 
            dc = Math.Log(br - bl);
            result = res24 * dc;
            abserr = Math.Abs((res24 - res12) * dc);
            res12 = 0.0;
            res24 = 0.0;
            for (i = 0; i < 13; i++)
            {
                res12 += (cheb12[i] * rg[i]);
                res24 += (cheb24[i] * rg[i]);
            }
            for (i = 0; i < 13; i++)
            {
                res24 += (cheb24[i] * rg[i]);
            }
            _130:
            result = (result + res24) * factor;
            abserr = (abserr + Math.Abs(res24 - res12)) * factor;
            goto _270;

            //  This part is executed only if b = br
            // *
            // *  Compute the Chebyshev series expansion of the following function:
            // *
            // *  f2 = (0.5 *(b+bl-a-a)+0.5*(b-bl)*x)^alfa *
            // *      f(0.5*(b-bl)*x+0.5*(b+bl))
            // 
            _140:
            hlgth = 0.5 * (b - bl);
            centr = 0.5 * (br + bl);
            fix = centr - a;
            fval[0] = 0.5 * f.SolveAt(hlgth + centr) * Math.Pow(fix + hlgth, alfa);
            fval[12] = f.SolveAt(centr) * Math.Pow(fix, alfa);
            fval[24] = 0.5 * f.SolveAt(centr - hlgth) * Math.Pow(fix - hlgth, alfa);
            for (i = 1; i < 12; i++)
            {
                u = hlgth * x[i - 1];
                isym = 24 - i;
                fval[i] = f.SolveAt(u + centr) * Math.Pow(fix + u, alfa);
                fval[isym] = f.SolveAt(centr - u) * Math.Pow(fix - u, alfa);
            }
            factor = Math.Pow(hlgth, beta + 1.0);
            result = 0.0;
            abserr = 0.0;
            res12 = 0.0;
            res24 = 0.0;
            if ((wgtfunc == 2) || (wgtfunc == 4))
            {
                goto _200;
            }

            //  wgtfunc = 1  (or 3)  
            DqchebClass.Dqcheb(x, fval, cheb12, cheb24);
            for (i = 0; i < 13; i++)
            {
                res12 += (cheb12[i] * rj[i]);
                res24 += (cheb24[i] * rj[i]);
            }
            for (i = 13; i < 25; i++)
            {
                res24 += (cheb24[i] * rj[i]);
            }
            if (wgtfunc == 1)
            {
                goto _260;
            }

            //  wgtfunc = 3  
            dc = Math.Log(br - bl);
            result = res24 * dc;
            abserr = Math.Abs((res24 - res12) * dc);
            res12 = 0.0;
            res24 = 0.0;
            for (i = 0; i < 13; i++)
            {
                res12 += (cheb12[i] * rh[i]);
                res24 += (cheb24[i] * rh[i]);
            }
            for (i = 13; i < 25; i++)
            {
                res24 += (cheb24[i] * rh[i]);
            }
            goto _260;

            //  Compute the Chebyshev series expansion of the following function:
            // *
            // *      f3 = f2 * log(0.5*(b-bl)*x+0.5*(b+bl-a-a))
            // 
            _200:
            fval[0] *= Math.Log(hlgth + fix);
            fval[12] *= Math.Log(fix);
            fval[24] *= Math.Log(fix - hlgth);
            for (i = 1; i < 12; i++)
            {
                u = hlgth * x[i - 1];
                isym = 24 - i;
                fval[i] *= Math.Log(u + fix);
                fval[isym] *= Math.Log(fix - u);
            }
            DqchebClass.Dqcheb(x, fval, cheb12, cheb24);

            //  wgtfunc = 2  (or 4)  
            for (i = 0; i < 13; i++)
            {
                res12 += (cheb12[i] * rj[i]);
                res24 += (cheb24[i] * rj[i]);
            }
            for (i = 13; i < 25; i++)
            {
                res24 += (cheb24[i] * rj[i]);
            }
            if (wgtfunc == 2)
            {
                goto _260;
            }
            dc = Math.Log(br - bl);
            result = res24 * dc;
            abserr = Math.Abs((res24 - res12) * dc);
            res12 = 0.0;
            res24 = 0.0;

            //  wgtfunc == 4  
            for (i = 0; i < 13; i++)
            {
                res12 += (cheb12[i] * rh[i]);
                res24 += (cheb24[i] * rh[i]);
            }
            for (i = 13; i < 25; i++)
            {
                res24 += (cheb24[i] * rh[i]);
            }
            _260:
            result = (result + res24) * factor;
            abserr = (abserr + Math.Abs(res24 - res12)) * factor;
            _270:
            return result;
        }
    }
}